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ABSTRACT 

We present the SpineWeb framework for the topological analysis of the Cosmic Web and the iden- 
tification of its walls, filaments and cluster nodes. Based on the watershed segmentation of the cosmic 
density field, the SpineWeb method invokes the local adjacency properties of the boundaries between 
the watershed basins to trace the critical points in the density field and the separatrices defined by 
them. The separatrices are classified into walls and the spine, the network of filaments and nodes in 
the matter distribution. Testing the method with a heuristic Voronoi model yields outstanding results. 
Following the discussion of the test results, we apply the SpineWeb method to a set of cosmological 
N-body simulations. The latter illustrates the potential for studying the structure and dynamics of 
the Cosmic Web. 

Subject headings: Cosmology: theory - large-scale structure of Universe - Methods: numerical - 
Surveys 



1. INTRODUCTION 

The large scale distribution of matter revealed by 
galaxy surveys features a complex network of inter- 
connected filamentary galaxy associations. This net- 
work, which has become known as the Cosmic Web 
(|Bond et al.lll996l ). contains structures from a few mega- 
parsecs up to tens and even hundreds of Megaparsecs of 
size. The weblike spatial arrangement of galaxies and 
mass into elongated filaments, sheetlike walls and dense 
compact clusters, the existence of large near-empty void 
regions and the hierarchical nature of this mass distribu- 
tion - marked by substructure over a wide range of scales 
and densities - are its three major characteristics. Its ap- 
pearance has been most dramatically illustrated by the 
recently produced maps of the nearby cosmos, the 2dF- 
GRS, the SPSS and the 2MASS redsh ift surveys (e.g. 
iColless et ai~l[200l Euchra et al. Il200l . 

In this paper we introduce the SpineWeb formalism 
for analyzing the structure and topology of the Cosmic 
Web. It identifies the sheets and filaments in the Cos- 
mic Web, along with the large underdense void regions, 
and their mutual connection into the Spine of the cosmic 
matter distribution. The method is bas ed on the Wa- 
tershed Transform fWST. iBeucher llTTlM? ; . and is largely 
free of user-specific parameters and artificial smoothing 
scale(s). Its output will enable the study of the physical 
properties and dynamics of the individual morphological 
components, along with their topology and hierarchical 
characteristics. 



1.1. the Cosmic Web 

The Cosmic Web is the most salient manifestation 
of the anisotropic nature of gravitational collapse, the 
motor behind the f ormation of structure in the cos- 
mos ([Peebles I ll980T ). N-body computer simulations 
have profusely illustrated how a primordial field of 
tiny Gaussian density perturbations transforms into a 
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pronounced and intricate filigree of filamentary fea- 
tures, dented by dense compact clumps at the nodes 
of the networ k (iColberg. Krughoff fe Connolly! [20051: 
iSprineel et al~l 120051) . The filaments connect into the 
cluster nodes and act as the transport channels along 
which matter flows into the clusters. 

Fundamental understanding of anisotropic collapse on 
cosmologic a l sca les came with the seminal study by 
IZeldovichl (|1970f ) . who recognized the key role of the 
large sca le tidal for ce field in shaping the Cosmic Web 
(also see llcke II 19731) . The collapse of a primordial cloud 
(dark) matter passes through successive stages, first as- 
suming a flattened sheetlike configuration as it collapses 
along its shortest axis. This is followed by a rapid evo- 
lution towards an elongated filament as the intermediate 
axis collapses and, if collapse continues along the longest 
axis, may ultimately produce a dense, compact and viri- 
alized cluster or halo. The hierarchical setting of these 
processes, occurring simultaneously over a wide range of 
scales and modulated by the expansion of the Universe, 
complicates the picture considerably. Recent state-of- 
the-art compjrtCT_e2£perhnents like the Millennium simu- 
lation (|Springel et al. I [20051 ) clearly show the hierarchi- 
cal nature in which not only the clust ers build up but also 
the fil amentary network itself (see lAragon-Calvo et al.l 
I2007H ). 

The Cosmic Web theory of iBond et all $1996) suc- 
ceeded in synthesizing all relevant aspects into a 
coherent dynamical and evolutionary framework. It 
is based on the realization that the outline of the 
cosmic web may already be recognized in the primor- 
dial density field. The statistics of the primordial 
tidal field explains why the large scale universe looks 
predominantly filamentary and why in overdense re- 
gions sheetlike membra nes are only marginal features 
(jPogosvan et al. 1 [1998). Of key importance is the 
observation that the rare high peaks, which will even- 
tually emerge as clusters, are the dominant agents for 
generating the large scale tidal force field: it is the 
clusters whi ch weav e the cosmic tapestry of filaments 
(IBond et al.l 119961: Ivan de Weygaert fc Bertschinger I 
119961: Ivan de Weygaert fc Bond I l2008al ). Thev cement 
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the structural relations between the components of 
the Cosmic Web and themselves form the junctions at 
which filaments tie up. This relates the strength and 
prominence of the filamentary bridges to the proximity, 
mass, shape and mutual orientation of the generating 
cluster peaks: the strongest bridges are those between 
the richest clusters that stand closely together and point 
into each other's direction. 

The emerging picture is one of a primordially and hier- 
archically defined network whose weblike topology is im- 
printed over a wide spectrum of scales. Weblike patterns 
on ever larger scales get to dominate the density field as 
cosmic evolution proceeds, and as small scale structures 
merge into larger ones. Within the gradually emptying 
void regions, however, the topological outline of the early 
weblike patterns remains largely visible. 

1.2. Closing in on the Cosmic Web 

Despite a large variety of attempts, as yet no gener- 
ally accepted descriptive framework has emerged for the 
objective and quantitative analysis of the geometry and 
topology of the Cosmic Web. The great complexity of 
both the individual structures as well as their connectiv- 
ity, the lack of structural symmetries, its intrinsic mul- 
tiscale nature and the wide range of densities that one 
finds in the cosmic matter distribution has prevented the 
use of simple and straightforward techniques. 

Historically, the quantitative analysis of the Cosmic 
Web has been dominated by a description in terms of 
statistical measures of clustering of galaxies and matter. 
While correlation functions have been the mainstay of 
the cosmological analysis of large scale structure, a di- 
rect interpretation in terms of the patterns and texture 
of the Cosmic web has largely remained elusive. Over 
the years a variety of heuristic measures have been for- 
warded to analyze specific aspects of the spatial patterns 
in the large scale Universe, but only in recent years there 
have been attempts towards developing complete descrip- 
tors of the intricate spatial patterns that define the Cos- 
mic Web. Nearly without exception these methods bor- 
row extensively from other branches of science such as 
image processing, mathematical morphology, computa- 
tional geometry and medical imaging. 

Noteworthy examples includ e filament detection with 
the help of the Candy model (jStoica et al. I 12005!) and 
wave let analysis of the Cosmic Web ([Martinez et al. I 
2005). Several methods seek to relate morphological 
features to singularities in the density field, usually in- 
voking information on the gradient and Hessian of the 
density field, or of the tidal field (s ee e.g.lSousbie et al. I 
2008aHAragon-Calvo et al.ll2007alfbl : IHahn et al Il2007al lbl: 
Bond et all 120091 ) . A classification scheme on the ba- 
sis of the manifolds in the tidal field - involving all 
morphological features in the c osmic matte r distrib u- 
tion - has been pre s ented by IHahn et all (|2007al rbl): 
iForero-Romero et al. I (|2008l) . However, its success may 
depend strongly on the correct choice of the smoothing 
scale. Another concept addressing the gradient and Hes- 
sian of the density field is that of the skeleton analy- 
sis, a direct appli cation of Morse theory to cosmological 
densi ty fields (see IColombi et al. 1120001 : iPogosvan et al. I 
2009). The skeleton formalism has been developed for 
the morphological analysis of the Megaparsec Cosmic 
Web, in redshift surveys like SDSS as well as in N- 



body si mulati ons ([Novikov et al.l l2006t ISousbie et al. I 
2008a. b, 2009). Its present implementation refers to 
features identified at one single specific scale and suf- 
fers from Gaussian smoothing. The multiscale nature of 
the cosmic matter distribution is explicitly addressed by 
the Multiscale Morphology Filter, which is based on a 
scale-space analysis of the Hessian of the density field 
lAragon-Calvo et ail (|2007al fbh to identify cluster, fila- 
ments and sheets on the scale where they are locally most 
prominent. 

1.3. Watershed and Cosmic Spine 

One technique that implicitly addresses the topology 
of the Cosmi c Web is the Waters hed Void Finder (WVF) 
developed bv lPlaten etall (|2007| ). The WVF is an appli- 
cation of the Watershed Transform for the identification 
of underdense basins in the megaparsec-scale matter dis- 
tribution. The Watershed transforms segments the den- 
sity field into isolated basins and delineates the bound- 
aries of cosmological voids. 

The method presented in this paper is the natural ex- 
tension of the Watershed Void Finder. It includes the 
watershed transform into a wider context as a frame- 
work for studying both the morphology and topology of 
the cosmic web and its various constituents. The result 
is the SpineWeb method, a complete framework for the 
identification of voids, walls and filaments. Via the prac- 
tical role of the watershed transform in computing the 
Morse complex it is intimately related to Morse theory, 
in which it finds its mathematical foundation. It is im- 
portant to note that although Morse theory can be used 
to describe the same topology traced by the SpineWeb 
method there is not a univocal relation between the two 
as the SpineWeb is based on the observed properties of 
the Cosmic Web. 

An important aspect of our method is that it is an 
intrinsically scale-free method, starting from a scale- 
free reconstruction of the density field. We use th e 
DTFE method of ISchaap fc van de Wevgaertl (|2000l ). 
which guarantees an optimal and unbiased representa- 
tion of the hierarchical nat ure and anisotropic morphol- 
ogy o f cosmic structure (see Ivan de Wevgaert k, Schaap I 
I2009L for an extensive description). Having guaranteed 
the capability of invoking a full scale- free Scale-space rep- 
resentation of cosmic structure, our watershed procedure 
not only traces the outline of filaments and sheets, but 
may also be extended towards doing so over a range of 
scales in order to address their hierarchical structure. 

1.4. Outline 

The principal rationale behind the SpineWeb analysis 
of the cosmic matter distribution is the interest in relat- 
ing the geometry of the matter and galaxy distribution 
in a more meaningful fashion to the underlying dynami- 
cal evolution. One particular aspect of this dynamically 
motivated disentanglement of structure is the attempt to 
identify various evolutionary stages of the tidally induced 
anisotropic collapse of structure in the Universe. 

In this paper we will focus specifically on the descrip- 
tion of the basic SpineWeb formalism, confined to a den- 
sity field sampled on a regular grid. We start by dis- 
cussing the topological background of our study in sec- 
tion^ focusing on the watershed transform, its connec- 
tion to the general context of Morse theory and the re- 
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Fig. 1. — Top left: a slice of 100 h~ 1 Mpc of side showing the density field computed from a N-body simulation. Top right: the same 
slice as the left panel but here showing the density field as a shaded landscape where high density regions correspond to ridges while 
underdense regions correspond to valleys. The light source used to shade the surface is located at the north-east. Bottom left: the 2D 
watershed transform computed from the 2D density field using a discrete-intensity bucket algorithm (black contours). Each individual 
valley is randomly colored for visualization purposes. Bottom right: the same landscape as the top right panel but here we also show the 
watershed transform superimposed as thick black lines. The watershed contour has been slightly smoothed for visualization purposes. 



lated issues of practical interest to the Spine Web formal- 
ism. The overall cosmological background of the struc- 
ture, formation and dynamics of the Cosmic Web fol- 
lows in section amongst others to establish the link 
between the structures identified by the Spine Web tech- 
nique and the filamentary i dentity of tidal b r idges in the 
theor y of the Cosmic Web fiZcldovich~1 ll970t iBond et al.l 
1996). The technical aspects of the Spine Web formalism 
are outlined in some detail in section Subsequently, 
the formalism is tested by applying it to two different 



classes of spatial particle distributions. The first testbcd 
concerns two simple heuristic Voronoi clustering models 
which model aspects of cellular and/or weblike spatial 
distributions. Visual and quantitative tests are described 
in section [5] The operation of Spine Web in a more re- 
alistic setting of a ACDM simulation is the subject of 
section [6l In this section we also stress the fundamental 
differences between a structural selection based on den- 
sity thresholds or one based on topological criteria. To 
illustrate the potential for analyzing cosmological struc- 
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tures, in section [7] we shortly describe three quantitative 
measures for the matter distribution in the ACDM simu- 
lation we used for testing. Finally, section [3] summarizes 
our results, and discusses prospects and further develop- 
ments of the Spine Web formalism. 



2. WATERSHED SEGMENTATION OF THE COSMIC WEB 

When studying the topological and morphological 
structure of the cosmic matter distribution in the Cosmic 
Web, it is convenient to draw the analogy with a land- 
scape (see fig. [TJ top row). Valleys represent the large 
underdense voids that define the cells of the Cosmic Web. 
Their boundaries are sheets and ridges, defining the net- 
work of walls, filaments and clusters that defines the Cos- 
mic Web (cf. top panels fig.[TJ. 



2.1. the Watershed Transform 

The watershed transform (WST) is one of the most 
common methods used in Image analysis for segmenting 
images into distinct patches and features. It is a concept 
defined within the context of mathematical morphol- 
ogy, a nd was first introduced by iBeucher fc Lantueioul I 
(1979). The basic idea behind the WST stems from geo- 
physics, where it is used to delineate the boundaries 
of separate domains, i.e. basins into which yields of 
e.g. rainfall will collect. The watershed transform is 
formed by the ridges and sheets surrounding the water- 
shed basins and includes a subset of all the critical points 
in the density field. 

The word watershed finds its origin in the analogy of 
the procedure with that of a landscape being flooded by 
a rising level of water. Suppose we have a surface in 
the shape of a landscape (cf. top right panel, fig. [T]). 
The surface is pierced at the location of each of the min- 
ima. As the water-level rises a growing fraction of the 
landscape will be flooded by the water in the expanding 
basins. Ultimately basins will meet at the ridges defined 
by saddle-points and maxima in the density field. The 
final result of the completely immersed landscape is a 
division of the landscape into individual cells, separated 
by ridge dams (see left bottom panel fig. [1} . 

2.2. A watershed search for voids 

The watershed transform was first introduced in a cos- 
mological context as an objective technique to identify 
and outlin e voids in the cosmic matter and galaxy dis- 
tribution (|PlatiniraD EQQ2 EliED H00l). Following 
the density field-landscape analogy, the Watershed Void 
Finder (WVF) method identifies the underdense void 
patches in the cosmic matter distribution with the water- 
shed basins. The method is parameter free in case there 
is no noise in the data. 

A major advantage of the WVF method is its inde- 
pend ence of assumptions on the shape and size of voids 
(see (jColberg et al. If2008h for a comparison of its perfor- 
mance with a variety of void finding algorithms). Sharing 
this virtue with a similar tessellat ion-based void finding 
method, ZOBOV (jNevrinck 112001 ) . WVF is particularly 
suited for the analysis of the hierarchical void distribu- 
tion expected in the commonly accepted cosmological 
scenarios. 



2.3. Watersheds and Landscape Gradients 

Extrapolating its application to other areas of interest, 
the implementation of the watershed transform may also 
be seen as a practical instrument for the segmentation 
of surfaces and volumes on the basis of the topological 
structure of the "landscape" /(x). To trace the topolog- 
ical structure of a field /(x), we need to investigate the 
structure of the gradient field of the landscape, V/ (for 
an exce llent introduction to computational topology, we 
refer to lEdelsbrunner I (|2010ft b 

2.3.1. Gradient Field and Integral Lines 

The gradient delineates a smooth vector field, which 
vanishes at critical points, 



V/(x k ) =0, 



(1) 



The integral lines or slope lines represent the flow along 
the gradient field V/ between the critical points. On the 
basis of these connections one may infer a variety of spa- 
tial segmenta t ions (see e.g.lCavlev 1118 59; Ma xwell 1I187C ; 
Eberlv 1119941 iFurst 112003: lEdelsbrunner et al.l l2003allb: 
Danovaro et al.l 120031 : iGvulassv et al. I I2005D . One par- 
ticular segmentation is the watershed transform, which 
segments the landscape / into regions of uniform local 
gradient behavior: the watershed basin j consists of the 
collection of points x that are closer in topographic dis- 
tance 7~(x, Yj ) , 



T(x,y,) = inf y |V/(7(a))|ds. 



(2) 



to the defining minimum y^ of the basin than to any of 
the other minima. In this definition the integral is the 
pathlength along the integral line, the line along whose 
path the tangent at each point is parallel to the local 
gradient V/. The watershed itself then consists of the 
ridge lines that delimit the boundaries between basins in 
the terrain. 

An illustration of the close link between the gradient 
field and structural features in the Universe is offered 
by the righthand panel of fig. [5J The image shows that 
the integral lines that define the boundaries of adjacent 
valleys are in fact the watersheds. It also reveals the inti- 
mate relationship between the critical points in the flow 
field and the nodes, filaments and voids in the landscape: 
maxima are found at nodes of the weblike network of 
watershed ridges, minima at the centers of the void cells, 
while saddle points are to be found at key locations along 
the ridges. Following this view, we see that the water- 
shed lines are the set of slope lines emanating from saddle 
points and connecting to a local maximum or minimum. 
Within this framework, saddle points have the crucial 
function of defining the sheets and filaments in the den- 
sity field through their connection to the maxima via the 
integral lines. 

Note that because the image in fig. [2] is a slice through 
a three-dimensional field, the identification between the 
structural elements and the critical points in the image 
is not entirely unequivocal (see below). Nonetheless, 
the principal observation is that the resulting weblike 
segmentation of space, and the corresponding boundary 
manifolds, contain the full information on its topological 
structure marked by sheets, filaments and nodes. 



Fig. 2. — Left: A slice of the density field shown as a shaded landscape with the watershed lines superimposed as black lines. Right: The 
zoomed area in the blue square of the left panel showing the slope lines (white lines) superimposed on the density field (gray background). 
The contour of the watershed transform is delineated by the thick black lines. 



2.4. Morse theory 

The vast majority of applications of the watershed 
transform concern the interior of the segmented regions. 
However, it is straightforward to extend its focus to other 
morphological components of the Cosmic Web, towards 
the delineation of the network of overdense ridges and 
walls which form the boundary manifolds of the cosmic 
density landscape. 

This can be directly appreciated by noting the close 
relation between the definition of the watershed trans- 
form and the more formal concept of the Morse com- 
plex. Morse theory is the mathematical framework for 
the analysis of the topological structure of manifolds, by 
relating it to smooth, C 2 -differentiable, functions defined 
on those spaces. Central to Morse theory are the loca- 
tion and nature of the critical points - minima, maxima 
and saddle points - and their mutual connection via the 
gradient-based integral lines. These determine the mor- 
phological features of the functional surface. 

Even tho ugh there are s ome differences between the 
two (see e.g. lGvuias sv 2008), the close similarity between 
the definition of the watershed transform and the con- 
cepts of Morse theory indicates that the computation 
of the watershed transform may be used as an efficient 
means of computing the various structural elements in 
a landscape dissected alon g the lines of Morse theory 
(M orsc l ll9^[MihTo71[l965l) . 

In a cosmological context, the skeleton f ormal ism 
(|Novikov et al.1 [200l ISousbie et al"~l l2008allbl . f2009h is 
also based on Morse theoretical concepts, via the gra- 
dient and/or Hessian of the density field. The approach 
followed in our Spine Web procedure involves the specific 
application of the watershed transform for the analy- 
sis and description of the topology of the Megaparsec 
Universe, following our introduction of th e concept in 
the context of cosmic density field analysis (jPlaten et al.l 



l2007f l. Implici tly, this results i n a pseudo- Morse se gmen- 
tation (see e.g. Gyulassy 2008; Edclsbrunner 2010), with 
the advantage of opening the path towards a fully hier- 
archical form alism. This intimate r elationship was also 
recognized bv ISousbie et al. I (|2009h . who used a proba- 
bilistic extension of the watershed technique in the latest 
version of the skeleton formalism. 

2.5. The Discrete Watershed Transform 

The implementation of the watershed transform in a 
large variety of scientific applications has to address a 
few important practical issues. A typical characteristic 
of most scientific images is their discrete nature. 

The discreteness concerns two aspects: the spatial dis- 
creteness, ie. the discrete number of intervals at which 
the image has been sampled (pixels/ voxels), and the dis- 
crete intensity levels at which the image has been sam- 
pled. 

Image discreteness creates a few complications for an 
accurate calculation of the watershed transform. It ren- 
ders it difficult to identify the existence and exact loca- 
tion of saddle points on the basis of a discretized local 
neighborhood. For the same reason, it is difficult to ac- 
curately extract slope lines. 

Several methods for the extraction of critical points 
have been developed in an attempt to alleviate the 
limitations imposed by the discreteness of images. 
Among these, the dis crete watershed transform algo- 
rithm (|Beucher I [1982) represents a simple and ele- 
gant formalism for identifying the watershed separatri- 
ces and can be shown to co nverge to the continuous case 
(|Naiman fc Schmitt 111994} ) . The procedure emulates the 
flooding of valleys or catchment basins in a (discrete) im- 
age representing a landscape. The points where two or 
more lakes converge are marked, and the algorithm con- 
tinues until all the pixels in the image have been flooded. 
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At the end of the process the image will be segmented 
into individual regions sharing a local minima, with the 
points that were marked as the dividing boundaries be- 
tween two or more valleys defining the watershed trans- 
form. 

A major asset of the intensity discretization is that it 
helps to remove faint features, and therefore also removes 
artefacts without the need of pre- or post-processing. 
Perhaps the greatest advantage is that discrete images 
allow the use of highly efficient algorithms but in general 
their use is limited to image segmentation since they give 
incorrect topologies. 

In the case of images with continuous (floating point) 
values one retains the option of computing the watershed 
transform directly from the continuous intensity image, 
in addition to the option of discretizising the intensity. 
On the basis of the continuous image, the watershed 
transform would delineate the topology more accurately 
than would be feasible on the basis of the the discrete- 
level representation. However, it would involve a sub- 
stantial increase in computational cost and of complexity 
of the code. 

3. SPINE OF THE COSMIC WEB 

The analogy between the watershed transform defining 
the boundary between underdense basins and the topol- 
ogy of the cosmic matter distribution is in itself one of the 
major justifications of the SpineWeb method presented 
in this study. Basic is the connection between the ele- 
ments that form the Spine of the Cosmic Web: walls, 
filaments, clusters and voids. 

• The Cosmic Web is an interconnected system of 
dense compact clusters, elongated filaments and 
tenuous sheetlike walls. Visible through the galax- 
ies, gas and dark matter populating these struc- 
tural features, the Cosmic Web theory (|Bond et al.l 
|1996[ ) teaches us that its topological outline was al- 
ready present in the primordial perturbation field 
out of which all structure arose. 

• All of the elements of the Cosmic Web are inter- 
connected. This is a crucial observation, which 
can be most readily appreciated b y studying high 
resolu tion N-body simulations (e.g. lSpringel et al. I 
I2005t ). Otherwise seemingly isolated objects usu- 
ally turn out to be connected to less massive struc- 
tures which become visible when assessing the mass 
distribution at a higher mass resolution. A tanta- 
lizing idea is that the galaxies found at the cen- 
ter of voids lie at the inters ection of tenuous intra- 
void dark- matter filaments. (Zitrin & Bro sch IfeOOijl : 
iPark fc Leell2009HStanonik et al. 112010ft . 

• Filaments are suspended between clusters or, de- 
pendent on scale, massive halo clumps. Their 
prominence and density may vary substantially, de- 
pendent on the mass, distance and alignment of the 
generating dark matter halos. However, the sheer 
presence of two matter clumps is already sufficient 
for the corresponding tidal force field to guarantee 
the topological presence of a filamentary bridge. 
Tenuous membranes permeate the space between 
adjacent filaments, and are part of the large wall 



which defines the boundary between two under- 
dense voids. The wall boundary is outlined by var- 
ious filaments, connecting each other at the cluster 
nodes. 

Following these observations, the Cosmic Spine is de- 
fined as the topological network of nodes, filaments and 
sheets along which the cosmic matter distribution on 
large Megaparsec scales has assembled (see fig-HJ). 

4. THE SPINEWEB PROCEDURE 

The key aspect of the SpineWeb procedure is that it 
exploits the intrinsic topological information contained in 
the Watershed Transform to delineate the Cosmic Spine. 
For the computation of the Cosmic Spine by means of 
the Watershed Transform it is necessary to address a 
few issues of practical importance. 

4.1. Density field 

In our current implementations of the SpineWeb pro- 
cedure, we apply the DTFE method to reconstruct the 
density field from the spatial particle distribution. 

The DTFE procedure produces a self- 
adaptive volume-filling density field on the ba- 
sis of th e Delaunay tessellation of the point 
distri bution (|Schaap fc van de Wevgaert 1 120001 : ISchaap I 
I2007T ). DTFE density (and velocity) fields have been 
found to optimally trace a hierarchical matter distri- 
bution at any resolution level represented by the point 
sample, while at the same time resolving the local 
anisotropies in the matter distribution. This high level 
of sensitivity to the topology of the matter distribution, 
makes DTFE ideally suited for the SpineWeb procedure 
(fo r an extensive description of the DTFE procedure, 
see Ivan de Wevgaert fc Schaap Il2009h . 

Given a spatial distribution of points, DTFE is based 
on the assumption that the density at the position of each 
point is proportional to the inverse of the total volume of 
the adjacent Delaunay tetrahedra, ic. to the volume of its 
contiguous Voronoi cell. Subsequently, the density field 
values at any location throughout the sample volume is 
determined by means of linear interpolation within the 
Delaunay tetrahedra of the corresponding Delaunay tes- 
sellation. Because a singular density determination at 
the central location of a voxel tends to introduce aliasing 
artefacts at high densities, we follow a slightly elaborate 
procedure. The density at each voxel of the "image" grid 
is determined on the basis of the DTFE density field sam- 
pled on a subgrid with a three times higher resolution 
and the density at each gridpoint set equal to the mean 
of the DTFE values at the 27 subgrid locations within 
the corresponding voxel. 

In the applications described in this study, we compute 
the density field values at the voxels of a regular cubic 
grid. We use a fast and efficient implementation of the 
DTFE algorithm based on the publicly available CGAL 
library 3 . 

4.2. Watershed Implementation 

The discrete watershed transform code we use in the 
SpineWeb procedure is an adaptation of the immersion 
and the topographical distance algorithms for floating 

3 www.cgal.org 
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point intensity values (see iRoerdink fc Meiister I l200Ct 
for a review). The code assumes a density map which 
is sampled on a regular grid. Our C code 4 computes the 
watershed transform from a double-precision 512 3 grid in 
just a couple of minutes on a regular linux workstation. 

In a first step, the code starts by finding and labeling 
the local minima in the density map, by identifying the 
voxels with the lowest density value among all their 26 
neighbors. These local minima are the seeds of the void 
valleys to be identified by the watershed transform. 

In a second step, we follow the topographical distance 
algorithm in order to obtain a fast segmentation of the 
space into locally connected underdense regions. For 
each voxel we identify the voxel among its 26 neighbors 
which has the lowest density. The maximum gradient 
paths are traced by iteratively connecting the voxels to 
their lowest density neighbor until the path reaches a lo- 
cal minimum. Subsequently, we assign the label of the 
corresponding minimum to the path. 

In the third step we extract the watershed transform 
itself, ie. the boundaries between the void regions. The 
pixels in the watershed transform are identified by means 
of a local immersion algorithm. First, we identify all the 
voxels that lie at the boundaries between two or more 
regions. Subsequently, this subset of voxels is sorted 
in density and the standard immersion algorithm is ap- 
plied. By following this two-step procedure, we avoid 
what would be the most expensive component of the al- 
gorithm. Instead of having to sort the complete density 
field, the sorting evaluations are restricted to the points 
in the watershed boundaries, a minor fraction of the com- 
plete volume. 

In a final step, each of the pixels in the watershed 
boundary is assigned a morphological label, following its 
identification as void, wall or filament element, according 
to the criterion expressed in equation [3] and as illustrated 
in fig. El 

While the watershed transform provides us with a 
highly efficient means of segmenting space into topolog- 
ically well-defined elements, in general the watershed al- 
gorithm tends to overdo the segmentation, creating too 
many regions and identifying small noise in the image 
rather than real features fsee lEdelsbrunner~ll2010t ). We 
circumvent this circumstance by preprocessing the den- 
sity or distance field via a Gaussian smoothing of u = 2 
voxels. 

4.3. from Watershed to SpineWeb 

From the analogy between the Cosmic Web and the 
watershed transform one can define, on the basis of the 
discrete watershed transform of a cosmic density field, 
a set of unique criteria to identify voids, walls and fila- 
ments. 

The criteria are based on the properties of the local 
neighborhood of all the points that comprise the discrete 
watershed transform. Instead of computing at any given 
point the local eigenvalues of the Hessian of the density 
field, one may simply resort to the entirely equivalent 
evaluation of the identity of the surrounding 26 neighbor 
pixels (for the three-dimensional situation) . By counting 
the number Af voids of adjacent watershed basins (voids) 

4 The code will be publicly released in an upcoming article. In 
the meantime it can be provided upon request. 



amongst these, it is straightforward to discriminate be- 
tween voxels which belong to a void, a wall or a filament 
by means of the following set of rules: 



AC, 



= 1, void 
= 2, wall 
> 3, spine 

(filament + clusters) 



(3) 



Note that in the present implementation we do not dis- 
criminate between filament or cluster node, and instead 
consider them to be part of the same spinal structure. In 
a future implementation we will include a density max- 
imum criterion which would allow us to find the cluster 
nodes amongst the spine voxels. In addition we can easily 
identify regions in clusters by using one of the common 
halo finders available like friends of friends and isolate 
them from the identified filaments. 

Our Spineweb technique exploits a purely local crite- 
rion for identifying the morphological nature of boundary 
pixels in the density field: it utilizes the full geometric 
structure of the watershed transform to limit the evalu- 
ation to the direct neighbours of each p oint. This dif- 
fers fr om the implementation followed by ISousbie et al. I 
(2009) who included a probability propagation scheme in 
the watershed flooding procedure from where the differ- 
ent elements of their skeleton were determined by finding 
their intersections between the corresponding peak and 
void patches. 




wall 




filament 



Fig. 3. — Local neighborhood around a voxel (green) inside a 
wall (left) and filament (right). The blue voxels indicate walls and 
the red voxels filaments. The light gray cubes here represent voxels 
inside voids. The voxel inside a wall has two adjacent voids inside 
its neighborhood (A and B) while the voxel inside a filament has 
in this case three adjacent voids (A, B and C). 



The above criterion is a purely and solely a topologi- 
cal one. By definition, walls are the regions between two 
adjacent voids. Filaments are to be found at the inter- 
section of three watershed basins, at the intersection of 
3 walls (see fig. [3j. The success of these criteria can be 
appreciated from the 3-D surface maps in fig. [5] and the 
comparison between density and spine maps in fig. 1101 

4.4. Image Grid Representation 

While a regular grid facilitates the computation of 
the watershed transform and the subsequent topological 
identification of the various boundary pixels (see fig. [3]), 
its simplicity may also involve a few possibly artefacts. 

The first artefact relates to the discrete nature of the 
voxels in the density field. As a results, the filaments 
and walls have an artificial thickness - even if they would 



be infinitesimally thin - which makes the look pixelated 
or jagged. A particularly good illustration of this is the 
Spine obtained for the Voronoi clustering model in fig. [5] 
(panel c) where we explicitly render individual voxels as 
cubes. In the asymptotic limit of inhnitcsimally thin vox- 
els the discrete watershed will converge to the continuous 
case. 

The second artefact concerns the anisotropic nature of 
the local neighborhood of each voxel: Each voxel has 6 
neighboring voxels at a distance d — 1 (in voxel units) , 16 

neighbors at d = y2 and 8 neighbors at d = V3. it also 
involves an angular neighbor distribution deviating sub- 
stantially from angular isotropy. A possible alternative 
would be to limit the neighborhood evaluation to the 
6 most direct neighbors. However, the poor sampling 
might lead to a considerable risk of missing important 
topological information. For two-dimensional images the 
solution would be more straightforward. The use of a 
hexagonal grid would involve equal distance for all neigh- 
bor pairs and a perfectly uniform angular distribution. 
Unfortunately, an equivalent perfect grid for the three- 
dimensional situation does not exist. However, the use o f 
Centroid Voronoi Tessellations (CVT. lDu et al. I (|1999f )) 
would certainly help to alleviate the main artefacts. 

4.5. Galaxy Spine Assignment 

Physically, filaments and walls are not infinitesimally 
thin structures. To identify the particles or galaxies at- 
tached to them, we therefore need to the define a (natu- 
ral) thickness which encloses these objects. 

In the applications described below, we account for 
this by applying the dilation morphological operator to 
the voxels labeled as filament and wall. The process 
increases the thickness of filaments and walls by one 
voxel and this procedure can be performed iteratively 
to further increase the thickness. The dilation opera- 
tor was applied first to voxels labeled as wall and subse- 
quently to pixels labeled as filaments following the num- 
ber of degrees of freedom in the local variation of the 
density field, i.e. first wall s and subsequently filaments 
(|Aragon-Calvo et al.ll2007bl ). In our particular case a sin- 
gle iteration with a 3 x 3 x 3 kernel provides a good result 
without excessively fattening the structures. 

5. SPINEWEB TEST: 
VORONOI CLUSTERING MODELS 

To test and quantify in an objective way the identi- 
fication of walls and filaments with the Spine Web pro- 
cedure we apply it to a few realizations of Voronoi 
Clustering models of the large s cale matter distribution 
(Ivan de Wevgaert fc Icke 1 1 19891 : Ivan de Wevgaert 1 12Q02L 
2010). 

The Voronoi clustering models are heuristic models for 
cellular spatial patterns which use the geome tric (and 
convex) structure of the Voronoi tessellation (|Voronoi I 
Il908t IOkabell2000l) to emulate the cosmic matter distri- 
bution. They offer flexible templates for cellular patterns 
and are easy to tune towards a specific spatial cellular 
morphology. This makes them very suited for studying 
clustering properties of nontrivial geometric spatial pat- 
terns. Because the location, geometry and identity of the 
various spatial components in Voronoi models are known 
precisely, they are ideal as testbeds for the SpineWeb 
procedure. Unless otherwise specified, the seeds of the 



tessellation usually involve a set of Poisson distributed 
pints. 

The Voronoi models use the corresponding tessellation 
for defining the structural frame around which matter 
will gradually assemble during the formation and growth 
of cosmic structure. Points are distributed within this 
framework by assigning them to one of the four distinct 
structural components of a Voronoi tessellation: 

• Void: 

regions located in the interior of Voronoi cells. 

• Wall: 

regions within and around the Voronoi cell faces. 

• Filament: 

regions within and around the Voronoi cell edges. 

• Clusters: 

regions within and around the Voronoi cell vertices. 

What is usually described as a flattened "supercluster" 
consists of an assembly of various connecting walls in the 
Voronoi foam, while elongated "superclusters" of "fila- 
ments" usually include a few coupled edges. Vertices are 
the most outstanding structural elements, corresponding 
to the very dense compact nodes within the cosmic web 
where one finds the rich clusters of galaxies. 

Among a variety of possible Voronoi clustering real- 
izations, two distinct yet complementary classes of mod- 
els are the most frequently used ones, the structurally 
rigid Voronoi Element M odels and the evolving kin ematic 
Voronoi models (see e.g. Ivan de Wevgaertll201fl for an 
extensive description) . Here we use one Voronoi Element 
Model, a composite of all four distinct components, and 
one realization of a Voronoi kinematic model. 

In the case of the Voronoi Element model, the walls, 
edges and vertices are infinitely thin, yielding a pure 
geometric discrete realization of the underlying Voronoi 
tessellation. The Voronoi kinematic model represents a 
more realistic situation in which the various structures 
are assigned a finite width. We first analyze the topo- 
logically cleaner configuration of the Voronoi Element 
model, to assess the performance of SpineWeb under op- 
timal conditions. Subsequently, we investigate whether 
in how far it can sustain this performance under less op- 
timal circumstances. 

5.1. Voronoi Element model 

Based on a uniform distribution of M cell seeds in a 
(periodic) box of size L, we start with a uniform dis- 
tribution of N particles throughout a (periodic) box of 
size L. The particles are distributed within the tessel- 
lation by projecting them - with respect to the seed of 
the cell in which they are originally located - onto the 
walls, edges or vertices surrounding their Voronoi cell. 
As a result each of the walls, filaments and vertices have 
a different density, although the density remains uniform 
within each of the individual elements. 

The Voronoi Element models we used for our test con- 
sisted of a set of particles distributed in a box of size L 
with 10 6 particles from which 70% resides in the walls, 
25% in the filaments and 5% in the clusters. We chose 
a realization completely devoid of particles in the inte- 
rior of voids. This distribution of particles was chosen in 
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Fig. 4. — Left: slice of the distance field computed from the Voronoi seeds (background). The particles located inside filaments, walls and 
nodes are also shown. Right: Voronoi model with density field reconstruction. A slice across the intensity of the density field is shown in 
orange scale in the background. The black lines correspond to the watershed transform and the particles are indicated by small diamonds. 
Blue particles are misclassifications (in all the categories) black particles are correct classified particles. 



order to obtain a more uniform sampling over the three 
morphologies compared to a more realistic distribution. 

5.1.1. Voronoi Distance Field 

For the purpose of this test, instead of basing the 
SpineWeb procedure on the reconstructed (and noisy) 
density field we use the knowledge of the underlying tes- 
sellation to define a clean distance field. 

The main idea behind the SpineWeb method, the iden- 
tification of morphological structures on the basis of their 
topology, does not depend on practical details of the den- 
sity field determination from a dataset of observed galaxy 
locations or computer simulations. In this respect, it is 
important to realize that the validity of the SpineWeb 
procedure can be tested and assessed on the basis of 
the topological structure of any field that is topologi- 
cally equivalent to the density field of the Cosmic Web. 
This indeed is true for the distance field, and any generic 
field marked by a monotonic increasing value from a field 
minimum towards its watershed transform. 

In the particular implementation described in this sec- 
tion, the distance field is defined as the Euclidean dis- 
tance from each particle to its closest Voronoi seed. Re- 
gions close to the cell centers have low values while re- 
gions in the planes and edges of the cell have large values, 
with the value gradually increasing along the direction 
from the Voronoi cell centers towards the projected loca- 
tion on the walls. Within the walls, the highest density 
is reached at the surrounding edges, ultimately peaking 
at the vertices. In the resulting distance field, small cells 
correspond to low field values, while the larger cells yield 
higher field values, particularly near their boundaries. 
Following this definition, the distance field emulates the 
range of densities encountered in the Cosmic Web. The 
finite volumes of Voronoi cells in periodic tessellations 
assures a convergence of the distance field on the bound- 



aries of data volume. 

Various definitions of distance fields might be used, 
largely dictated by the specific questions at hand. An 
example of one such possibility would be to normalize the 
distance field on the basis of the distance to of the point 
to its second closest Voronoi nucleus, or the distance of 
its projected location on the corresponding Voronoi wall. 
The resulting distance field would reach a value of unity 
at the walls of the Voronoi cells. However, we opted to 
use the distance field with no normalization since it gives 
a better representation of the large dynamical range of 
densities encountered in the Cosmic Web. 

5.1.2. Distance Field Realization 

For the Voronoi Element test model, we determined 
the distance field on a cubic 512 3 grid. For each of the 
pixels in the grid we identify the closest nucleus, among 
the set of M generating nuclei. The pixel is assigned the 
value of its Euclidean distance. 

The resulting field is shown in fig. |4j It depicts the 
distance field itself, in a planar section through the 3- 
D box, by means of a greyscale map. The corresponding 
particles located in the filaments, walls and nodes, within 
a narrow strip around the sectional plane, are superim- 
posed on the image. 

5.2. Voronoi Kinematic Model 

A nearly equivalent Voronoi clustering model realiza- 
tion is a Voronoi Kinematic Model. Here initially ran- 
domly distributed particles move away from their expan- 
sion nucleus - ie. the closest nucleus in whose Voronoi 
cell they are located - by a universal expansion rate 
(van de Wevgaertll200i MM . see e.g.). 

When particles reach the wall shared between their 
expansion center and the second closest nucleus, their 
motion is constrained to their path within the wall. This 




10 






Fig. 5. — The SpineWeb method applied to a Voronoi distribution, a) original particles lying at the edges of the Voronoi cells (filaments), 
b) original particles lying at the faces of the Voronoi cells (walls), c) Pixels inside filaments (red) and walls (blue) identified with the 
SpineWeb method, d) recovered particles lying at the edges of the Voronoi cells (filaments), e) recovered particles lying at the faces of 
the Voronoi cells (walls), f) particles erroneously identified as particles in filaments. The box shown here contains 1/64 of the original box 
volume. 



continues until they reach one of the edges delimiting the 
wall, upon which they proceed along the edge towards 
their ultimate location, that of one of the vertices at the 
ends of the edge. 

The model simulation box has a length of 141h _1 Mpc, 
in which we find 180 Voronoi cells. In total, the box 
contains N = 2097152 particles. Originally distributed 
randomly throughout the box, we move them until 16.5% 
of the galaxies reside in the walls, 28.7% in the filaments 
and 51.3% at the cluster nodes. Unlike the Voronoi Ele- 
ment Model described in the previous section, the voids 
remain populated with a diluted random distribution of 
3.5% of all void galaxies. 

The different morphological structures in this Voronoi 
kinematic model are also assumed to have a finite phys- 
ical width. Particles within the walls, edges and vertices 
are assumed to have a Gaussian distribution perpendic- 
ular to these structures. The width for each of these 
morphologies is set to a = 1.0h _1 Mpc. Even though 
the topological properties of the pure Voronoi Element 
Model and this Voronoi kinematic model are practically 
equivalent, the finite width and more organic develop- 
ment of the Voronoi kinematic models represent a spatial 
density field which more closely emulates that encoun- 
tered in galaxy redshift surveys and in N-body simula- 
tions of structure formation. 



5.2.1. Voronoi Model Density field 

The SpineWeb performance test for the Voronoi kine- 
matic model is based on the density field of the particle 
distribution. We use the DTFE method to reconstruct 
the density field from the point distribution on a 256 3 
cubic grid within the simulation box. 

FigurelH righthand panel, shows a 2-D section through 
the reconstructed density field, which is represented by a 
color map. The low-density patches and high-density 
structures clearly outline the interior region of voids, 
while the high-density regions correspond to the walls 
and filaments in the overall spatial pattern. 

5.3. SpineWeb identification 

Following the determination of the distance field for 
the Voronoi Element Model, and the density field for 
the Voronoi kinematic model, the SpineWeb procedure 
proceeds to identify the Spine of the particle distribution. 
Following the identification, we assess the fraction of false 
and real spine detections in both Voronoi models. 

5.3.1. Detection Rates: definition 

Quantitatively, we assess the detection rate of the 
SpineWeb procedure by determining the ratio of real and 
false detections of filament and wall particles. 
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The detection rate i? roal is denned as the fraction of 
original wall or filament particles that are also identified 
as such by the Spine Web technique, 

(4) 

1 * original 

In this equation, AT rea] is the number of wall or fila- 
ment particles that have also been identified as such, 
and iV original is the total number of particles that in the 
Voronoi model physically belong to a wall or filament. 
Along the same lines, we measure the filament and wall 
contamination rate R {B , lso following the definition, 

N talM is the number of particles recognized as wall or 
filament particle, but in reality having a different mor- 
phological identity. 

The results for the the detection and contamination 
rates of walls and filaments in both the Voronoi Element 
model as well as the Voronoi kinematic model are listed 
in table Q] 



TABLE 1 

Recovered particles per morphology 



Structure 




-Rfalso 


Walls dis 


0.93 


0.15 


Spine dis 


0.91 


0.03 


Walls DEN 


0.76 


0.32 


Spine den 


0.87 


0.13 



TABLE 2 

Ratio of real and false recovery rates per morphology. 

The top half corresponds to the results from the 
distance field (dis), while the bottom half concerns the 
results for the DTFE reconstructed density field (den). 
For definition of i? HEAL and R false see eqn.[2]and \E\ 



5.3.2. Spine of the Voronoi Element Model 

In order to remove small-scale spurious variations, the 
distance field is smoothed with a Gaussian filter of a = 
2 voxels. Of the smoothed distance field we compute 
the watershed transform. Following this, the Spine is 
determined. 

The Spine Web results for the Voronoi Element Models 
are illustrated in figure [5j Visual inspection of the figure 
provides a good impression of the virtues and perfor- 
mance of the procedure. Comparison between panels (a) 
and (d) shows that the genuine filament particles in the 
Voronoi model (a) are identified with a convincing accu- 
racy by the Spine Web procedure (d). The same is true 
for the successful identification of the Voronoi wall parti- 
cles (panel b) and the Spine Web identified wall particles 
(panel e). Interestingly, the particles that Spine Web er- 
roneously identify as filament particles while in fact they 
are wall particles, shown in panel (f), clearly delineate 
the original filamentary web. It is most most likely a 
result of the discrete resolution of the distance field grid. 

The Spine Web reconstructed spine, shown in panel (c) 
of fig. [5l allows a clear and transparent assessment of the 



topological structure of the Voronoi web. It shows the 
identified wall voxels by means of blue blocks, while the 
filamentary voxels are indicated as red blocks. The walls 
form a continuous network of connecting surfaces, with 
filaments delineating the intersections between the walls. 

From table Q] we can see that 93% of particles in walls 
and 91% of particles in filaments are also identified as 
such, while the contamination rate of walls and filaments 
are 15% and 3%. The false identities can usually be 
ascribed to the jagged nature of the voxels (see panel (c), 
fig.®. 

The simple and idealized example of the Voronoi El- 
ement model shows the intrinsic potential and power 
of the SpineWeb procedure. Solely on the basis of the 
topology of the matter distribution, and independent of 
a density threshold or any other arbitrary parameter, it 
manages to determine its correct morphological segmen- 
tation. 

5.3.3. Spine of the Voronoi Kinematic Model 

The situation for the Voronoi Kinematic Model, 
marked by a more noisy particle distribution and a less 
idealized density field reconstruction, is less straightfor- 
ward yet more representative for realistic circumstances. 

Fig. [4] (righthand panel) shows, superimposed on the 
density field color map, the Voronoi particle distribu- 
tion as well as the watershed transform. The latter is 
shown by means of the black lines. Particles identified 
as wall or filament particles are indicated by small dia- 
monds. The watershed transform is able to identify most 
of the voids and their boundaries. The image shows that 
the SpineWeb procedure managed to closely follow the 
location of the edges and faces in the original Voronoi 
tessellation. 

Some minor deviations are detected at small scales, the 
result of Poisson noise in combination with the discrete 
nature of the sampled density field. By construction, 
our filaments and walls are only one voxel thick (in this 
Voronoi model realization this is equal to ~ l/i Mpc). 
As a result, we may miss the particles inside a given 
structure because of small-scale variations in the den- 
sity field translating into variations in the watershed 
transform. This effect can be appreciated from the 
miss-classified (blue) particles in fig. HJ mostly clustered 
around the boundaries between filaments and filaments. 

To determine the detection and contamination rate of 
the SpineWeb calculation, we compare the identity of 
the Voronoi wall and filament particles, which are a pri- 
ori known from the model generation, with that of the 
classification on the basis of the reconstructed watershed 
segmentation. To this end, we assign the particles within 
a radius of 2a ~ 2/i _1 Mpc from a wall or a filament in the 
spine (watershed) segmentation to that particular struc- 
ture: a particle is identified with a wall when it it lies 
within a 2a distance from two different watershed cells. 
Although the detection rate results are less forthcoming 
than for the pure Voronoi element models, they remain 
convincing (see lower half of table [U We find a detection 
rate of 76% for wall particles, as opposed to a contami- 
nation rate of ~ 32%. Filaments are better recognizable, 
which may be understood from the 87% detection rate 
of filament particles, as opposed to a mere 13% misclas- 
sification rate. 
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Fig. 6. — Volume rendering of a thick slice of the density field in a ACDM simulation in a box of 200h Mpc. Left: the full DTFE 
reconstructed density field. Right: the density field inside an isosurface at <5 = 1. 



6. SINGLE-SCALE ACDM SPINE 

To test the SpineWeb method in a more realistic and 
challenging situation we have applied it to a cosmolog- 
ical N-body simulation. It concerns a ACDM universe 
simulation inside a box of 200 h~ x Mpc, restricted to the 
dark matter particles. Initial conditions were generated 
on a 512 3 grid with Vt m = 0.3, Cl A = 0.7, ct 8 = 0.9 
and h = 0.73. For the pr imordial pertu r bation s, we 
use the transfer function of iBardeen et al.l ()1986l ) with 
a shape paramet er T = 0.21 following the definition of 
ISugivamal (|1995l) . After having set up the initial condi- 
tions, we follow the subsequent gravitational evolution to 
the present time using the public N-body code Gadget2 
(jSpringel et al. Il2005h . 

From the same initial conditions, we also generated 
additional lower-resolution versions of 256 3 , 128 3 and 
64 3 particles were genera ted, following the "a veraging" 
prescription described in (jKlypin et al. ll2~00lD . For the 
single-scale analysis in this paper, which focuses on the 
largest filaments and walls in the particle distribution, it 
is sufficient to analyze the low-resolution 64 3 dataset. 
The lower resolution corresponds to a cut-off scale of 
~ 3h~ 1 Mpc in the initial conditions, sufficient for the 
analysis focusing on the large scale structure. The higher 
resolution datasets are used for visualization purposes. 

6.1. Density field morphology 

From the final particle distribution we compute the 
density field on a cubic grid of 512 voxels per dimen- 
sion using the DTFE method (see sect. 14. 1|) . Figure |6] 
(lefthand panel) depicts a volume rendering of the den- 
sity field in a thick slice through the simulation box. It 
shows that DTFE manages to follow the intricacies of 
the weblike structures in great detail, over a range of 
scales: it reproduces the correct geometry of the various 
features. 

An interesting example of structural complexities in 



the displayed region is the cluster at the lefthand side of 
the slice. A full 3D visualization of the system shows that 
the filaments entering the cluster define several semi- 
planar structures, all sharing the cluster as their com- 
mon node. Lower isodensity contours reveal even more 
of the tenuous walls, even though at such low density 
levels we need to take into account that the image gets 
easily confused by spurious interloping features. 

A frequently used approach for delineating structural 
features such as filaments and walls is to assign a specific 
density range to each morphology. Filaments or walls are 
singled out by selecting the regions that have a density 
within the corresponding density range. In the righthand 
panel of fig. |6]we show the density field inside an isoden- 
sity surface at S = 1, for clarity superimposed on top 
of a white background. Because a density value 5 = 1 
is roughly comparable to ty pical values encountered i n 
filaments and walls (see e.g. lAragon-Calvo et al.ll2~007h . 
these contours roughly define the boundaries of the fila- 
ments and the clusters embedded, along with the tenuous 
walls suspended between the filaments. 

Although the isodensity surfaces provide good insight 
into the overall distribution of matter, one immediate 
observation is that the attempted pure density selec- 
tion of filaments is not very succes sful. As was pointed 
out by Arago i>Calvo et "all (|2007f) filaments and walls 
are characteri zed by a rather br oad range of densities 
(also see e.g. lHahn etaLl l2007al ). The broad cluster, 
filament, wall and field density ranges are also mutu- 
ally overlapping over a siz eable density range (also see 
lAragon-Calvo et al.|[2010cT ). 

An immediate repercussion of the overlapping density 
ranges is that is almost impossible to decide purely on 
the basis of a density criterion whether a certain location 
belongs to a cluster node, filament or wall. This may be 
directly appreciated from the truncated density map in 
fig. |U Although the image shows a substantial degree of 
filamentary structure, comparison with the full density 
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Fig. 7. — Three-dimensional watershed transform of the density field. The cube shows the pixels that compose the watershed transform, 
from which the Spine is extracted. Several slices cut across the simulation box show the watershed (white lines) delineating the density 
field (blue-green background). The three dimensional nature of the watershed network is evident. 



field shows that it discards the pattern of lower density 
filaments. Also, it does not manage to disentangle the 
highly concentrated agglomeration of filaments near the 
massive cluster at the central lefthand side of the box. 
Moreover, throughout the whole volume it is rather dif- 
ficult to see which locations would belong to a filament 
and which ones to a wall. 

6.2. Cosmic Spine and Cellular Morphology 



Following the computation of the density field, we com- 
pute its watershed transform. The resulting segmenta- 
tion of the density field into its watershed basins is illus- 
trated in figure [7] 

The watershed basins are to be identified with the void 
regions in the cosmic matter distribution. To get a better 
idea of its spatial structure and the connections between 
the various structural components, we slice through the 
watershed field at regular intervals along the a;-axis. This 
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yields a sample of yz-slices through the simulation box. 
Figure [7] shows a sequence of five consecutive yz-slices, 
with the watershed segmentation shown as white lines su- 
perimposed on top of the density field (blue-green level) 
map. It is straightforward to appreciate the correspon- 
dence between the watershed segmentation lines and the 
underlying density field. While the segmented cube in 
figure [7] emphasizes the characteristic cellular nature of 
the Cosmic Web, the yz slices reveal the close relation- 
ship between voids, walls and filaments. 

The 2-D yz slices show the strong correlation between 
the density field and the watershed transform. The wa- 
tershed lines trace the high-density ridges and regions 
in the density field, occasionally bridging their lower- 
density connecting parts. On a more global scale, we also 
notice that the higher density regions contain a higher 
number of distinct and smaller cells than the more mod- 
erate or underdense areas. This translates into a more 
complex local network of filaments and walls. The oppo- 
site effect occurs in underdense regions. These are mainly 
characterized by large symmetrical voids, surrounded by 
relatively simple wall-filament environments. 

Cosmic voids are immediately recognized as large 
empty cells in the watershed transform. The walls in 
the cosmic matter distribution are visible as the bound- 
aries between two adjacent watershed cells, while the fil- 
aments are found at the intersection of these walls. The 
considerable variety of sizes and shapes of voids is most 
readily visible in the pattern of watershed lines in the 
yz-slices. Even though we know that on stereological 
grounds, lower-dimensional sections tend to exaggerate 
the size distribution of the full 3-D distribution, the com- 
parison with the void basins in the 3-D box does confirm 
the impression of the diversity of the void population. It 
also underlines the significance of topological Spine Web 
analysis: the void distribution is a direct reflection of the 
complexity of the dyna mical processes which are forming 
and shaping the v oids ([Sheth fc van de Weveaertl 12004 
iPlaten et al1l2008t) . 

6.3. Cosmic Spine: Filaments and Walls 

The final step in the SpineWeb method is to identify 
and label the voxels that correspond to the spine (fila- 
ments and clusters) and walls, following the SpineWeb 
criteria specified in equation [3] 

An insightful impression of the intricacy of the full 
three-dimensional network of filaments and walls is pre- 
sented in fig. El The top two frames show the wall-like 
(blue) and filamentary (red) regions separately. Clearly 
outstanding is the percolating nature of the filamentary 
network and the complex of connecting sheets. 

The appearance of a uniform width, at places consid- 
erably in excess of the local width of the density field 
contours, is a result of our choice to show, for visualiza- 
tion purposes, a uniform smoothed outline. The plot- 
ted isosurface of both walls and filaments is obtained 
by filtering the mask defined by all wall voxels with a 
Gaussian kernel of a = 2 voxels. The Gaussian smooth- 
ing radius corresponds roughly to the average width of 
ps 2/i -1 Mpc of filaments and walls, as w e found in a pre- 
vious study (|Arag6n-Calvo et al.l l2007ft . One may get 
an impression of the corresponding variation in density 
and width along the spinal structures by inspecting the 
two bottom frames of fig. 1101 where we have superim- 



posed density contour levels onto the embedding spinal 
contours (left: walls; right: filaments). 

The top-left panels in figure |S]shows a network of com- 
plex sheets. Instead of a regular "planar" geometry, on 
small scales the walls have a curved appearance marked 
by an irregular surface. To a considerable extent this 
reflects their inhomogeneous internal mass distribution, 
itself a result of their hierarchical buildup. The irregular 
convoluted shapes are found on all scales, although the 
walls do have a slightly more regular semi-planar geom- 
etry on larger scales. 

Also the filamentary structures reflect their inhomoge- 
neous internal mass distribution, even though the applied 
smoothing has tended to diminish the contrast between 
e.g. massive clumps and tenuous moderate or lower den- 
sity parts of the filaments. As a result, over most of its 
outline the filamentary edges look semi-linear. Nonethe- 
less, occasionally we can recognize rather twisted con- 
figurations, and at numerous locations we can recog- 
nize the bulging presence of massive clusters. Taking 
the nodal junctions as the endpoints demarcating an in- 
dividual filament, a first analysis shows that short fila- 
ments tend to be more straight than longer ones. This 
is entirely in line with the trend predicted by the Cos- 
mic Web theory, and is in agreement with a correspond- 
ing analysis of N-body simulations (jBond et al.l 119961 : 
iColberg. Krughoff fc ConnoIbTl [20051 : ?). 

6.4. Morphological Connections 

The close mutual relationship between the walls and 
filaments is immediately clear when inspecting the su- 
perposition of the wall-like and filamentary web in the 
bottom (right) frame of figure [H With the filaments 
defining an interconnected web, the walls fill the spaces 
in between the filaments, together forming a "watertight" 
complex of membranes surrounding a system of voidlike 
cavities. The zoom-in onto one specific region centered 
around a particularly intricate branching of filaments em- 
anating from a node, in fig. [HI highlights the complex con- 
nections that may occur in the Cosmic Web. At least four 
filaments appear to originate from a core region at the 
confluence of five walls. In the branches we can clearly 
recognize the bulging imprint of massive clusters. The 
zoom-in also nicely shows the small-scale heterogeneity 
of the walls' surfaces. 

6.5. Cosmic Spine versus Density Selection 

To compare the morphological SpineWeb segmentation 
with the corresponding density field, in figure [10] we have 
depicted the matter distribution in a central slice through 
the simulation box. The two top panels are images of the 
density field (see discussion sect. I6.1| ). while the bottom 
panels present the density field within morphologically 
segregated regions. The two frames at the bottom show 
the density field inside semi-transparent surfaces enclos- 
ing the wall features (bottom left) and the filamentary 
features (bottom right). 

6.5.1. Stereological Considerations 

When assessing the bottom frame of fig. QUI we have 
to take into account that isolated features observed in 
these slices are the result of the finite thickness of the 
slice. Particularly noteworthy for the filaments in the 
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Fig. 8. — Surfaces enclosing the voxels which are identified as belonging to walls (blue, top left) and filaments (red, top right) within a 
cubic region of 50/i -1 Mpc. The bottom frame shows how in the same region both morphological components are connected and intertwined. 
The latter forms a nice illustration of the intimate relationship between filaments and walls. For visualization purposes the surfaces are 
smoothed with a Gaussian kernel of <r = 2 voxels. 




Fig. 9. — Zoom-in onto the cosmic spine in a subregion of the 50h 1 Mpc, highlighting the intricate connections between wall surfaces 
(blue), filamentary edges and cluster nodes (red). 

righthand frame, they are artefacts of the finite width of the depicted slice. Because the orientation of filaments 
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Fig. 10. — Filaments and walls identified with the SpineWeb algorithm. Top left: Volume rendering of the density field inside a subbox 
of the simulation. Top right: density field contained inside an isosurface at <5 = 1. Only the density inside the isocontour is plotted. 
Bottom left: density field and semi-transparent isosurfaccs delineating walls. Bottom right: density field and semi-transparent isosurfaccs 
delineating filaments. Both wall and filament mask have been smoothed for visualization purposes with a Gaussian filter of a = 2 pixels. 



in the cosmic spine with respect to the slice is random, 
their intersection with the slice will differ. Dependent on 
the intersection angle, it may vary from their full length - 
in case they are nearly entirely embedded within the slice 
- to a mere point in case they run perpendicular to the 
slice. The resulting impression is one of a semi-irregular 
distribution of shorter and longer "stubs" , which indeed 
we find back in the bottom righthand frame. 

The two-dimensional geometry of walls tends results in 
linear intersection with the depicted slice. This produces 
the fully percolating network of (intersection) edges seen 
in the bottom lefthand frame. Note that occasionally the 
orientation of a wall is so favorable that its intersection is 



not a one-dimensional edge but instead consists of a slab 
comprising a major fraction of the wall. In the most ex- 
treme circumstance, the wall is entirely embedded within 
the slice so that it remains visible in its entirety. 

6.5.2. Morphological Structure of Density Features 

While the first superficial impression might be that the 
isocontour map of the density field (see fig. fTUl top-right 
panel) is richer in detail and structure than the filamen- 
tary and sheetlike morphologies in the bottom frames, a 
few important observations need to be made. 

An important contrast between the isodensity contour 
maps and the filamentary and sheetlike networks defined 
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by the Spine Web procedure is that between the rather 
discontinuous nature of the (thresholded) density maps 
and the fully percolating spinal structure. This is most 
clearly visible when zooming in regions with large den- 
sity contrasts, of which filaments close to the infall region 
are a good example. The massive cluster complex visi- 
ble at the left of the box forms a nice illustration. The 
filamentary extensions connecting to the cluster are iden- 
tified by the Spine Web procedure (lefthand panel fig. fTU)) . 
along with the sheetlike membranes of which they form 
the boundary (righthand panel fig. HOp . It would be very 
challenging for traditional density-based filament detec- 
tion techniques to trace filaments near cluster-filament 
interfaces. The density in the infall regions of clus- 
ters tends to increase dramatically, rendering a density- 
based criteri on to determine the local m orphology rather 
cumbersome (|Arag6n-Calvo et al.ll2007t) . 

7. ANALYSIS WALL & FILAMENT SAMPLE 

In this section we present a few quantitative mea- 
sures of the voids, walls and filaments extracted by the 
Spine Web technique. The results concern the single-scale 
analysis of our 64 3 particle ACDM simulation described 
in the previous section (see sec. IB]). 

7.1. Density distribution 

The density field of the 64 3 simulation was computed 
on a regular 512 3 grid, no smoothing was applied to the 
filament and wall masks. Figure [TT] shows the density 
field distribution computed for the complete simulation 
box, voids, filaments and walls. The distribution of den- 
sities can be roughly described as log-normal with the 
main difference between morphological environments be- 
ing the position of the peak of the distribution. Numer- 
ous studies have shown that a gravitationally evolving 
matter distribution, starting from Gaussian initial condi- 
tions, tends to attain a lognormal density distribution to- 
wards more advanced quasi-linear stages ([Coles fc Jones I 
U99lHNevrinck et al.ll2009HPlaten |[2009h . Our results in- 
dicate that this remains true for each of the individual 
morphologies. 

Voids have the lowest densities followed by walls and 
filaments respectively. It is important to note that the 
network of filaments found by our method contains also 
the clusters which act as the nodes of the wall-filament 
network. This affects the right tail of the distribution and 
the computed moments. The median is a much better 
estimator than the mean in all cases given the effect of the 
clustering of matter into halos in our sampling schema. 

Table [3] shows basic statistics of the density field char- 
acterized by the Spine. We present the statistic com- 
puted at two different grid resolutions (256 3 and 512 3 
voxels) as a simple convergence test. The numbers in 
table [3] show that the mean and median densities of the 
different morphologies are largely similar for the two dif- 
ferent resolutions. This is quite different for the vol- 
ume occupancy, in particular for the walls, filaments and 
clusters. This is a direct reflection of the resolution- 
dependent finite thickness assigned to filaments and 
walls, i.e. the voxel size. The volume fraction of voids is 
less sensitive to the grid resolution. We may understand 
this in terms of the Spine Web invariants: void occupancy 
scales by volume, wall occupancy by surface area and fil- 
aments by length. The latter is therefore most sensitive 
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Fig. 11. — Normalized density distribution for all voxels in the 
simulation box (thick gray line), void voxels (solid line), wall voxels 
(dotted line) and filament voxels (dashed line). 

to density field resolution, voids least. 



TABLE 3 

VOIDS, FILAMENTS AND WALLS: DENSITIES AND VOLUMES 



Structure 





(1 + 5)256 




-^m,256 




mean 


median 




voids 


0.81 


0.56 


0.825 


walls 


1.68 


0.94 


0.158 


spine 


3.46 


1.58 


0.015 




(1 + 5)512 




•7"m,512 




mean 


median 




voids 


0.86 


0.58 


0.889 


walls 


1.85 


0.93 


0.103 


spine 


4.47 


1.57 


0.006 



TABLE 4 

Basic statistics of the density field in voids, walls and 
filaments, for density grids of 256 3 (top) and 512 3 
(bottom) voxels. quoted are the mean and median of the 
density field, (1 + <5), in the individual morphologies, and 
the volume fraction tv,m of each of them. 



An important observation is also the considerable over- 
lap between the pdf's of the different morphologies. 
While distributions peak at clearly different places, there 
is a large overlap between all density distributions. This 
degeneracy in the density distributions between mor- 
phologies explains why a given isodensity contours does 
not manage to isolate one specific morphology, but will 
invariably include regions belonging to other morpholo- 
gies too. 

For our purpose, most significant are the differences 
between density distributions inside filaments and walls. 
Perhaps most remarkable is the sizeable overlap between 
densities in the void fields and those in filaments. The 
density distribution inside voids is almost identical to the 
overall density distribution. This is not surprising given 



18 




m ,„-i 



10" 



walls 



filaments 



_i i 1111 



_ i i i i '''' 



_i i i i ■ ■ ■ 



10 



100 



h" 1 Mpc 



Fig. 12. — Box-counting dimension of the filamentary spine 
(black dashed line) and walls (black solid line) of the Cosmic 
Web. For comparison we show the curves for ideal one and two- 
dimensional objects (gray dashed and gray solid lines respectively). 

the fact that even though voids are extremely underdense 
they occupy most of the space in the Cosmic Web. The 
difference between both distributions occurs at the high 
density tails where the clusters lie. 

7.2. Minkowski- Bouligand dimension 

We performed a preliminary scaling analysis of the 
identified filamentary and wall-like networks in the ana- 
lyzed ACDM N-body simulation. To this end, we have 
determined for each of these networks the Minkowski- 
Bouligand dimension D MB - or box counting dimension 
- formally defined as: 



lim 



logjV(e) 
o log(l/e) 



(6) 



In this expression, we count the number N(e) of boxes 
of (infinitesimal) size e required to fill or cover the set 
of points belonging to the filamentary or wall-like web. 
In practice, we divide the simulation box into subboxes 
of size s and count the number N(s) of subboxes that 
contain at least one voxel labeled as filament or wall. By 
repeating this evaluation for several box sizes, and deter- 
mining the scaling index N(s) oc s~~ D , we obtain an esti- 
mate of the Minkowski-Bouligand dimension. One may 
visualize this by plotting the count N(s) versus the size 
s if the boxes, preferentially in a logarithmic diagram. If 
indeed characterized by a single fractal dimension, the 
resulting curve would be characterized by one slope. In 
practice, the structural patterns tend to be more com- 
plex, manifesting itself in scaling curves that cannot be 
characterized by a single uniform slope. 

Figure [T2] shows the Minkowski-Bouligand dimension 
computed for the wall and filament networks. We also 
show two (grey) lines with slope -1 and -2 as a reference 
indicating the cases of a pure one and two-dimensional 
objects. The slope of the curves for filaments and walls 



differ considerably at small scales. Filaments behave like 
one dimensional lines up to scales of 3 — 4 h~ 1 Mpc after 
which point the absolute magnitude of the slope of the 
curve increases from -1 to -3 at scales of approximately 
10 /i _1 Mpc. In the case of the wall network we see a sim- 
ilar behavior with walls having a clear two-dimensional 
nature at scales smaller than 3 — 4 h~ 1 Mpc. The tran- 
sition point in the curve of figure [12] provides a good 
indication of the scale at which filaments and walls start 
joining each other forming an interconnected network. 
At this point their dimension is no longer 1 (filaments) 
or 2 (walls) but a higher value reflecting the complex- 
ity of the network of filaments and walls that form the 
Cosmic Web. 

7.3. LSS complexity and local density 

Another measure of the local complexity of the network 
of filaments and walls is presented in figure [T3] where we 
show the mean number of cells labeled as filament or wall 
inside boxes of 8 /i~ 1 Mpc size as a function of the mean 
density inside the boxes. 

Low values indicate very simple local configurations 
while large values reflect complex environments. At first 
glance this may seem straightforward as increasing ex- 
cursions sets of the density field have a similar behav- 
ior. However, the filaments and walls we identify are 
one- voxel thick so their voxel count correlates with their 
length and surface area respectively. For a given fixed 
volume larger counts indicate more intricate filament and 
wall systems. 

We find a trend between the density and the complex- 
ity of the environment. Highly dense boxes tend to con- 
tain more structures than underdense boxes. The regions 
in the vicinity of massive clusters are a good example of 
complex neighborhoods defined in a locally overdense re- 
gions while the large voids define relatively simple wall 
and filament structures. The spread about the clearly 
visible mean trend is quite substantial. One of the main 
reasons is the restriction in the scale of the analysis, 
which leads to a confusion of intrinsic scales and counting 
of fainter structures together with more significant ones. 
A proper multiscale analysis, the subject of our following 
paper, will take this into account. 

8. CONCLUSION AND FUTURE WORK 

The Spine of the Cosmic Web is the cosmic web's 
framework, consisting of the network of filaments and 
walls and their connections at the clusters nodes. In this 
study we present a topological technique based on the 
discrete watershed transform of the cosmic density field 
for the identification and characterization of voids, walls 
and filaments. Our method is closely related to a vari- 
ety of concepts from computational topology, and has a 
strong mathematical foundation in Morse theory of sin- 
gularities. 

The Spine Web method is ideally suited for morpholog- 
ical and dynamical studies of the Large Scale Structure. 
Amongst others, it will allow a better insight into the 
formation and dynamics of the anisotropic filamentary 
and wall-like structures in the Large Scale Universe. An- 
other immediate application is in addressing the question 
whether and which influences the large scale environment 
has on the halos and galaxies that are forming within 
their realm. 
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Fig. 13. — Number of cells labeled as wall (left panel) or filament (right panel) inside boxes of 8 h x Mpc of side as a function of the 
mean overdensity inside the 8 h~ 1 Mpc box. The number of cells is normalized with the mean count of all the 8 h _1 Mpc boxes. 



As a first test of its viability, we applied our method 
on a set of heuristic Voronoi clustering models. The 
Spine Web procedure succeeds in reconstructing the orig- 
inal properties of the cellular galaxy distribution. In the 
implementation presented in this work, we effectively re- 
strict ourselves to a single spatial scale determined by the 
voxel scale of the regular grid on which the density field 
is sampled. In a forthcoming paper we will discuss the 
effect of the multiscale nature of the matter distribution. 
The scale-space formulation of the Spine Web method will 
enable us to identify fainter features in the density field 
and establish their connections with other objects into 
a truly hierarchical weblike pattern. In other words, it 
provides an effective way towards characterizing the hi- 
erarchy of structures in the Cosmic Web. 

A crucial aspect of the watershed transform and of 
our method is the definition of local neighborhood. In 



the case of regular grids the immediate neighborhood 
of 26 pixels is arguably the best option. However, for 
unmeshed data such as galaxy surveys and N-body 
simulations, other neighborhood definitions offer a 
better choice. Among these, the Voronoi contiguous 
cell defined by the Delaunay tessellation of the point 
distribution represents a promising option. In the third 
paper of this series we will present the result of a 
Delaunay implementation of the Spine method. 
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